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Abstract 

The paper shows an application of the Smoothed Particle 
Hydrodynamics (SPH) for the numerical modeling of 
engineering problems involving rapid evolution over time as 
large strains and gradients, heterogeneity, deformable 
contours, mobile material interfaces and free surfaces. 

Following a Lagrangian approach, the continuum is 
discretized by means of a finite number of material particles 
carrying physical properties and moving according to 
Newton's equations of the classical physics. Spatial 
derivatives of a variable at a point are approximated by 
using the information on the neighboring particles based on 
the kernel approximation. 

The paper provides a brief description of the basics of the 
method along with some numerical aspects concerning 
boundary treatment and time integration scheme. 
Furthermore, some more details are provided about the 
recent improvements achieved with the aim of performing 
future SPH simulations involving underwater explosion for 
bottom sediment resuspension and flushing in an artificial 
reservoir. 

Numerical examples are illustrated and discussed 
concerning the early 2D test cases investigating the basic 
features of gas explosion; and obtained results show 
thateven if some improvements are required to overcome 
the model limitations, the SPH method is promising to 
reproduce the impulsive dynamics of the underwater 
sediments. Besides that, a feasibility study has been 
scheduled in order to investigate possible applications in 
innovative fields as the Enhanced Geothermal Systems for 
geothermal energy production. 
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Introduction 

The core idea of a meshfree method is to obtain a 



discretization of the continuum through a set of 
arbitrarily distributed nodes (or particles) that are lack 
of a connective mesh and can be adapted to possible 
topological and geometrical changes. 

With respect to traditional grid-based approaches, a 
meshless method allows tracing the deformation 
undergone by the material without excessive 
degradation of numerical results (owing to conflicts 
between mesh and physical compatibility) and high 
computational effort (e.g. adaptive mesh refinement). 

When the nodes assume a physical meaning (i.e. they 
represent material particles carrying physical 
properties such as mass, momentum etc.) the method 
is said a meshfree particle method and follows, in 
general, a Lagrangian approach. 

Among the different meshfree particle methods, the 
Smoothed Particle Hydrodynamics (SPH) was 
originally developed as a probabilistic model to 
simulate astrophysical problems [Lucy, Gingold & 
Monaghan], and later modified as a deterministic 
meshfree particle method and applied to continuum 
solid and fluid mechanics [Monaghan 1992, Monaghan 
1992a] because the kinematics and dynamics of the 
liquid particles, responding to Newton's equations of 
the classical physics, could be described in analogy 
with the simulation of the collective movement of 
astrophysical particles at large scale. 

According to standard SPH, a continuous physical 
quantity A(x), defined on the spatial domain Q as a 
function of the position vector x, and its spatial 
derivatives at the z-th material point, is approximated 
by using the information on the neighboring particles 
based on the kernel estimate. 

This procedure adopts a kernel function W(r, h), which 
is continuous, non-zero and depends on the modulus 
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of the relative position r= \xi- Xj I of the neighboring j- 
th particle falling within a circular space (spherical in 
3D problems) with radius 2/z, where h is generally 
referred to as the smoothing length (Fig.l). 



Radius of the compact 
support K-h 



Kernel function 




Particle /-th 



Particle /-th 



domain of 
particle /-th 



FIG. 1 TYPICAL REPRESENTATION OF PARTICLE 
DISCRETIZATION AND KERNEL FUNCTION 

The SPH approximation of the field function A(x) 
originates from the concept of integral representation: 

A(x) = \j(x-x')A(x')dn x , (1) 

In Eq.l, the Dirac delta function 8 is replaced by the 
kernel function leading to the kernel approximation: 

(A(x)) = \w(r,h)A{x')dn x , 

Q 

Considering the set of material particles representing 
the discretized continuum, the discrete form of the 
Eq.2 is obtained through the so called particle 
approximation that leads to the estimate of the field 
function A at the i-th particle: 

(A(X;)) = X— A(xJ%,*) (3) 

j=i Pj 

The summation in Eq.3 is extended over the N- 
neighboring particles, having volume AVj = rrijlp, 
falling within the compact support (or influence 
domain) Qi of the z-th particle. 

In a similar fashion, it can be demonstrated that the 
particle approximation of the function derivative can 
be obtained by shifting the differential operator on the 
kernel; and two alternative expressions are commonly 
adopted in fluid mechanics [Li & Liu, 2004]: 

(V ■ A(x,. )) = — f> . [A(x . )- A{x t )]• V W( rij , h) 
Pi j=t 



(V-A(x,.)) = A .£m. 



Vw(r ip h) 



(4) 



Pj Pt 

By applying the SPH interpolation, the Lagrangian 



form of the Navier-Stokes equations for a weakly 
compressible viscous fluid can be transformed into a 
system of ordinary differential equations that, by 
adopting the equations (4) and replacing W(ny, h) with 
Wij, are written as: 
'Dp\ 
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The additional term IT/ in Eq.5 is the so called 
Monaghan artificial viscosity [Monaghan 1992a] often 
introduced instead of the physical viscosity for 
numeric stability: 
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There are several advantages that can be obtained 
from such an approach: representation of the 
evolution of both free-surfaces, moving-interfaces and 
breaking waves, becomes more simple to face 
[Manenti 2008]; treatment of large deformation and 
shock problems becomes a relatively easier task [Di 
Monaco]; the particle tracking along with the relevant 
field variables can be obtained by numerical solution 
of the discretized set of governing equations in 
Lagrangian form [Manenti 2011]. 

Numerical Aspects 

The solution strategy of a meshfree particle method 
follows a pattern similar to a grid-based method. 

The computational domain is divided into a finite 
number of particles, followed by the numerical 
discretization of the system of partial differential 
equations according to the procedure described in the 
previous section and the resulting ordinary 
differential equations are solved through any stable 
time-stepping algorithm [Monaghan 2005]: here a first 
order explicit numerical scheme is used and a cubic 
spline function is adopted for kernel representation 
[Di Monaco]. 

The obtained velocity field allows one to update the 
particle position x and to compute the density field by 
means of the continuity equation (1); the pressure pi at 
each point is generally calculated through the equation 
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of state for a weakly compressible fluid and then 
smoothed out: 



: Pot + -(Pi-P t) 



av. w. 



(7) 



Solid boundaries are treated by means of the semi- 
analytic technique [Di Monaco]. Each portion of the 
solid contour, contributing to the mass and 
momentum equations of the generic z-th particle, is 
replaced by a fluid region extending beyond the 
boundary and treated as a material continuum with 
uniform velocity (m = ui), and hydrostatic pressure 
distribution (Fig.2). 

A typical term for boundary contribution in the 
balance equations is: 

CJC 2 J n ($,$)da 

with: Q =/(/>„ u„u fc/ v) (8) 

C 2 =f(Vr) J n (W) = Q 3 ^r n dr 

In Eq.8 a=f(3, cp) denotes the solid angle under which 
the z-th particle sees the portion of the solid boundary 
intersected by its sphere of influence and the integrals 
J n (n=l, 2, 3) depends on the boundary's geometry and 
can be computed analytically. 



Fluid-dynamic field 




extended fluid region 

FIG. 2 SKETCH OF SOLID BOUNDARY TREATMENT (2D CASE). 

Modeling Gas Explosion 

The explosion process of a high explosive (HE) 
material is characterized by a violent oxidation 
involving a chemical compound and an oxidizer; since 
the internal energy of the products is lower than the 
one of the reactants, a large amount of heat (say 



reaction heat) is quickly released [Cooper]. 

Even if such a phenomenon develops at very high 
speed of reaction, in the early phase it is characterized 
by two distinct inhomogeneous zones: a detonation- 
produced explosive gas and a non-oxidized explosive, 
between which a very thin layer exists representing 
the front of a reacting shock wave (detonation wave) 
advancing with a characteristic velocity U. 

Anyway, in several applications, the detonation speed 
can be assumed indefinitely high and the HE charge 
completely transformed into gaseous products whose 
expansion can be analyzed by considering the Euler 
equation for an inviscid fluid and assuming adiabatic 
process [Manenti 2008, Morris]. 

As a result, the viscous contribution at the right hand 
side of the linear momentum Eq.5 is neglected, and a 
balance equation for the gas internal energy is 
introduced: 



De 
Dt 



1 N 



Pi 



Pi Pi 



(u ( -n,)-VW 8 



(14) 



The state equation given by eq.15 is used for adiabatic 
transformation. 

P i ={r-l)p i e i < =y{y-l)e = - (15) 

Pi 

The kernel function adopted in subsequent analyses is 
a cubic spline, while the time integration is carried out 
through an explicit numerical scheme derived from 
the symplectic algorithm [Manenti 2010]. 

Numerical Examples 

This section illustrates the numerical results 
concerning two basic 2D SPH simulations of gas 
explosion: in the first of which a circular charge 
surrounded by a water crown expands inside a 
squared box with rigid walls; while in the other a 
squared charge expands inside an underwater 
sediment layer. The numerical model previously 
described is adopted. 

Gas Explosion 

In the following, the numerical results are shown 
concerning the expansion process of a HE detonated 
gas and the circular charge in vacuum is firstly 
considered (see Fig. 3); subsequently the charge is 
surrounded by a water crown (see Fig. 4) while in both 
cases the process is confined within a rigid box. As 
previously reported, the modulus of the detonation 
velocity U is assumed to be indefinitely high with 
respect to the gas kinematic, thus the explosive charge 
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is assumed completely detonated. 



FIG. 3 EXPANSION OF DETONATED GAS IN VACUUM INSIDE 
A RIGID BOX; AXES SCALE IN METERS 

Table 1 summarizes the relevant model parameters 
adopted in subsequent computations. 

When considering the vacuum gas expansion, at the 
initial time, 20 particles are placed in the radial 
direction while 60 particles are positioned along the 
tangential direction at the external boundary, resulting 
in a total number of 1200. The particle position, 
velocity (modulus) and pressure are depicted in Fig. 3 
at time intervals of 10 us and axes labels are in meters. 

The charge centre represents a singular point since no 



particle is placed on it and this explains the lower 
pressure. 

The expansion process reflects theoretical expectation 
until t = 40 us: past that time some irregularities in the 
pressure distribution at the outer boundary of the 
gaseous mass appear; therefore such a fact should be 
connected with the lack of information owing to the 
small number of neighbors in the interaction domain 
of external particles. 

table 1 Principal model parameters adopted for gas explosion 
computations 



MODEL PARAMETERS 


lo 


interp articles distance 


0.005 


m 


7] 


smoothing length 


1.3 ho 


m 


Po 


water ref . density 


1000 


kg/m 3 


P Y o 


gas ref. density 


1630 


kg/m 3 


sO 


spec, detonation energy 


4.29E+06 


J/kg 


Xa 


speed of sound 


5.0E+04 


m/s 


(Xm 


artificial viscosity param. 


0.2 




Pm 


artificial viscosity param. 


10.0 




s m 


velocity smoothing 


0.2 




Yn 


water state equation param. 


1.4 - 7.0 




7r 


gas state equation param. 


1.4 





Such non-physical behavior is however avoided if a 
surrounding medium is considered for confinement of 
the explosive charge, as illustrated in the following. 

Fig. 4 shows the expansion of a circular charge 
surrounded by a water crown and confined in a rigid 
squared box with length of 0.30 m and the simulation 
is carried out considering the same value of the state 
equation parameter for both gas and water (i.e. 
ycl yw=l, see Tab. 1). 

The upper left-hand panel displays the initial 
configuration and the position of the gauges for 
pressure detection on the transversal and diagonal 
directions and continuous green line denotes the rigid 
box walls. The central and the right-hand upper panels 
show particles position, velocity modulus and 
pressure at time t = 0.07 ms when the water impacts 
against the box walls. 

The lower panel shows pressure distribution at initial 
time (f=0.0ms) and at the impact time (f=0.07ms): in the 
latter a pressure wave is reflected by the box wall and 
propagates backward along the transversal direction 
with a pressure peak of about 1 GPa. 

The simulation ends at 0.2ms: after the gas and water 
particles have completely expand occupying the whole 
box internal volume. In addition, symmetrical jets 
originate from both transversal and diagonal 
directions thus pumping the gas toward the box center 
and producing a contraction of its volume (not 
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displayed in Fig. 4). gamma water constant in the state equation (i.e. yc / yw 

Fie. 5 shows the results obtained when the value of the * S mcrease< ^ to ^-0- 
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FIG. 4 EXPANSION OF DETONATED GAS SURROUNDED BY A WATER CROWN INTO A RIGID BOX; AXES SCALE IN METERS (yclyw = 1) 
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FIG. 5 EXPANSION OF DETONATED GAS SURROUNDED BY A WATER CROWN INTO A RIGID BOX; AXES SCALE IN METERS (yGlyw=Q2) 
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The particles dynamics is rather similar to that one 
described in the case yG / yw = V. anyway as now the 
water compressibility modulus is greater than the gas 
and this produces reflected pressure waves with 
higher peaks and celerity.The pulsation frequency of 
the gas expansion-contraction cycle is increased with 
respect to the case with yc / yw=l. 

Underwater Explosion Inside non-Cohesive Sediment 

In this numerical test, a squared HE charge is placed 
below the surface of a sediment layer at the bottom of 
a water tank, which represents a first attempt toward 
the investigation of a multiphase shock problem that 
becomes important in the development of an 
alternative methodology for siltation control in an 
artificial reservoir by means of combined use of 
explosive charges and a flushing manoeuvre (Manenti, 
2012). 

The principal parameters for the considered materials 
are summarized in Table 2. The adopted value of the 
sound speed has been calculated through the last 
equation in (15) assuming the gamma constant yG and 
the initial specific internal energy eo of the HE 
explosive. The compressibility modulus of each 
medium has been determined by assuming Cs=155 m/s 
and the density of the corresponding material. During 
the expansion, the water-sediment mixture has been 
simulated as first approximation, like a viscous fluid 
and this assumption appears reasonable based on the 
rapid increase of the pore water pressure in the solid 
matrix around the explosive charge, as discussed in 
the following. 



TABLE 2 MODEL PARAMETERS FOR WATER-GAS EXPLOSION 
COMPUTATIONS 



MODEL PARAMETERS 




Water 


Sediment 


HE explosive 


P [kg/m 3 ] 


1000.0 


1750.0 


1630.0 


ju [Pa s] 


0.0 


0.0 


0.0 




2.4 


2.4 


1.4 


eo [J/kg] 


252.0 


145.0 


4.29e+06 


h [m] 


0.06 


0.06 


0.06 


c s [m/s] 


155.0 


155.0 


155.0 


a M [-] 


2.0 


2.2 


1.0 


6>[-] 


0.2 


0.2 


0.2 



Figure 6 shows the obtained results. The upper panel 
displays the initial configuration and the distribution 
of the materials. The frame at t=260 us shows the shock 
wave propagation after the HE gas expansion, 
additionally, the initial value of eo=4.29e+05 J/kg for 
the specific internal energy has been assumed. 



FIG. 6 UNDERWATER EXPLOSION IN A NON-COHESIVE 
SEDIMENT DEPOSIT 
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The frame at t = 720 us shows the formation of a 
typical pinnacle just above the charge, and the 
pressure inside the sediment and water layers falls to 
values around 0.5 MPa and great part of the initial 
internal energy is transformed into kinetic energy, as 
illustrated in the lower panel showing the distribution 
of the velocity modulus and direction. 

Experimental investigations should be carried out in 
order to provide a data base to validate the SPH code; 
anyway the obtained results seem to be quite 
acceptable from a qualitative point of view and be in 
agreement with theoretical expectations. 

Note that there are several zones of the domain that 
are characterized by a coarse particle resolution, 
especially in the rarefied gaseous mass, which could 
cause the reduction of neighbours to very low values 
thus lowering the accuracy of the numerical 
computation. 

A particle refinement technique, such as that one 
proposed in [Liu] could be helpful for solving the 
above problem. 

Application Perspectives in Geophysics 

The advanced SPH applications previously described 
have been focused on underwater explosion modelling 
since this agreed with the target of the reference 
developing project. However, the good results 
obtained encourage the investigation of other possible 
fields of application of this innovative technique. One 
of the most promising application seems to be the 
study of the artificial fractures induced in warm dry 
rocks in order to create advantageous conditions for 
geothermal energy production by means of the so 
called Enhanced Geothermal System [Moia]. In this 
frame, a preliminary feasibility study has been 
scheduled with the aim of investigating the effects of 
pulses generated by high pressure fluid jets or 
microexplosions in a geological formation, also 
coupling a local SPH approach with a generalized 
CFD problem discretization. 

Conclusions 

An advanced application of the Smoothed Particle 
Hydrodynamics method for the numerical modeling 
of underwater explosion problems has been illustrated 
in this paper. 

The basic features of the numerical model adopted to 
simulate underwater expansion of a HE gas have been 
illustrated. 



The proposed results have shown that the physics of 
the investigated problems could be simulated with an 
adequate degree of accuracy for engineering 
applications; however, it requires some improvement 
to solve intrinsic criticisms, such as controlling the 
local particle refinement that could help to enhance 
the numerical accuracy in those regions where gas 
rarefaction takes place. 

In addition, further applications in geothermal energy 
production field are under investigation. 

ACKNOWLEDGMENT 

This work has been financed by the Research Fund for 
the Italian Electrical System under the Contract 
Agreement between RSE S.p.A. (previously ERSE) and 
the Ministry of Economic Development-General 
Directorate for Energy and Mining Resources 
stipulated on July 29, 2009 in compliance with the 
Decree of March 19, 2009 



List of Symbols 


A 


physical field function (scalar or vector) 


Cs 


speed of sound 


D/Dt 


material derivative 


e 


specific internal energy 


h 


smoothing length 


ho 


initial interparticle distance 


AV 


particle volume 


m 


mass 


n 


dimension of the physical space 


N 


neighboring particles 


V 


pressure 


po 


reference pressure 


Tij= 1 Xi-Xj 


1 modulus of the relative distance ve 


U 


detonation wave characteristic velocity 


Vr 


radial unit vector 


8 


gravitational acceleration vector 


u 


velocity vector 


Uij 


relative velocity vector 


Ub 


velocity vector of the solid boundary 


X 


position vector 


Xij 


relative position vector 


dV 


elementary volume 


OCM, Pm 


constants of Monaghan artificial viscosity 


8 


Dirac delta function 
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y state equation parameter 

s fluid compressibility modulus 

<p, 3, r spherical coordinates 

ju dynamic viscosity 

v cinematic viscosity 

II Monaghan artificial viscosity 

p density 

p G gas density 

po reference density 

9p pressure smoothing coefficient 

W kernel smoothing function 

dQ elementary volume of the continuum 

Q spatial domain 

/2/ compact support (or influence domain) of the 
z-th particle 
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